! ----------------------------
!
! Main program.
!
! ----------------------------
program resultsDisplay
  implicit none
  
  integer*8 m
  integer*8 numberOutputs
  character*(*) fileName
  
  parameter(fileName = '590000_600000_slaved_scheme.txt')
  parameter(m = 22)
  parameter(numberOutputs = 1000000)
  call displayMainLoop(fileName, numberOutputs, m)
end program resultsDisplay


! ----------------------------
!
! Main Loop.
!
! ----------------------------
subroutine displayMainLoop(fileName, numberOutputs, m)
  implicit none
  
  integer*8 m
  integer*8 option
  integer*8 numberInputs
  integer*8 numberOutputs
  character*(*) fileName
  complex*16 VELOCITY(m, numberOutputs)
  
  option = 1
  numberInputs = 0
  
  do while (option .ne. 0)
    print *, '1  Structure Functions'
    print *, '2  3th Order Moment vs 3th Order Structure Function'
    print *, '3  Shells Values'
    print *, '4  Velocity Map'
    print *, '5  Scaling Exponent'
    print *, '6  Energy'
    print *, '7  Velocity'
    print *, '8  Shells Modulus'
    print *, '9  Shells Correlations'
    print *, '0  Exit'
    print *
    print *, 'Option: '
    read *, option
    numberInputs = numberInputs + 1
    
    if (option .ne. 0) then
      call initPGPlot()
      
      if (numberInputs == 1) then
        call readVelocity(VELOCITY, fileName, numberOutputs, m)
      end if
    end if
    
    if (option == 1) then
      call displayStructureFunctions(VELOCITY, numberOutputs, m)
    else if (option == 2) then
      call display3OrderMomentVs3OrderStructureFucntion(VELOCITY, &
                                                        numberOutputs, m)
    else if (option == 3) then
      call displayShellValues(VELOCITY, numberOutputs, m)
    else if (option == 4) then
      call displayVelocityMap(VELOCITY, numberOutputs, m)
    else if (option == 5) then
      call displayScalingExponent(VELOCITY, numberOutputs, m)
    else if (option == 6) then
      call displayEnergy(VELOCITY, numberOutputs, m)
    else if (option == 7) then
      call displayShellsModulusThroughTime(VELOCITY, numberOutputs, m)
    else if (option == 8) then
      call displayEachShellModulus(VELOCITY, numberOutputs, m)
    else if (option == 9) then
      call shellsCorrelations(VELOCITY, numberOutputs, m)
    end if
    
    if (option .ne. 0) then
      call PGEND
    end if
  end do
  
  return
end subroutine displayMainLoop


! ----------------------------
!
! Display the structure functions of order p(from 2 to 10) with a log scale y
! and a linear scale in x
!
! ----------------------------
subroutine displayStructureFunctions(VELOCITY, numberOutputs, m)
  implicit none
  
  integer*8 m
  integer*8 order
  integer*8 numberOutputs
  real*4 X(m)
  real*4 Y(m) 
  real*8 STRUCUTURE_FUNCTION(m)
  complex*16 VELOCITY(m, numberOutputs)
  
  order = 10

  call structureFunction(STRUCUTURE_FUNCTION, VELOCITY, order, numberOutputs, m)
  call getLog10Values(STRUCUTURE_FUNCTION, m, X, Y)
  call setGraphViewAdjusted(X, Y, m, '\fin', '\fiStructure Functions', '', 20)
  call outLines(X, Y, m, 1)
  
  do order = 2, 9
    call structureFunction(STRUCUTURE_FUNCTION, VELOCITY, order, &
                           numberOutputs, m)
    call getLog10Values(STRUCUTURE_FUNCTION, m, X, Y)
    call outLines(X, Y, m, 1)
  end do
  
  return
end subroutine displayStructureFunctions


! ----------------------------
!
! Compares the Third-Order Moment with the Third-Order Structure Function.
! Log scale in y and linear in x.
!
! ----------------------------
subroutine display3OrderMomentVs3OrderStructureFucntion(VELOCITY, &
                                                        numberOutputs, m)
  implicit none
  
  integer*8 m
  integer*8 numberOutputs
  real*4 X(m) 
  real*4 Y(m)
  real*4 XLR(2)
  real*4 YLR(2)
  real*8 THIRD_ORDER_MOMENT(m)  
  real*8 STRUCUTURE_FUNCTION(m)
  complex*16 VELOCITY(m, numberOutputs)
  
  call thirdOrderMoment(THIRD_ORDER_MOMENT, VELOCITY, numberOutputs, m)
  call structureFunction(STRUCUTURE_FUNCTION, VELOCITY, 3, numberOutputs, m)
  call getLog10Values(STRUCUTURE_FUNCTION, m, X, Y)
  call setGraphViewAdjusted(X, Y, m, '\fin', &
                            '\fi< |u\dn\u|\u3\d >, \gP\dn\u/k\dn\u', '', 20)
  call outLines(X, Y, m, 1)
  call getLog10Values(THIRD_ORDER_MOMENT, m, X, Y)
  call outLines(X(1:m-1), Y(1:m-1), m-1, 2)
  call getLinearRegressionPoints(X, Y, 5, 17, m, XLR, YLR)
  call outLines(XLR, YLR, 2, 1)
  
  return
end subroutine display3OrderMomentVs3OrderStructureFucntion


! ----------------------------
!
! Displays the resulted zeta-p function with the corresponding scaling exponents
!
! ----------------------------
subroutine displayScalingExponent(VELOCITY, numberOutputs, m)
  implicit none
  
  real*4 X(m)
  real*4 Y(m)
  real*4 XR(9)
  real*4 YR(9)
  real*4 slope
  real*4 yIntercept
  real*8 STRUCUTURE_FUNCTION(m)
  integer*8 m
  integer*8 order
  integer*8 numberOutputs
  complex*16 VELOCITY(m, numberOutputs)
  
  do order = 2, 10
    call structureFunction(STRUCUTURE_FUNCTION, VELOCITY, order, &
                           numberOutputs, m)
    call getLog10Values(STRUCUTURE_FUNCTION, m, X, Y)
    call getLinearRegressionData(X, Y, 5, 17, m, slope, yIntercept)
    
    XR(order - 1) = REAL(order)
    YR(order - 1) = ABS(slope / REAL(LOG10(2.d0)))
  end do   
  call setGraphViewAdjusted(XR, YR, 9, '\fin', '\fi\gz\dp\u', '', 2)
  call outLines(XR, YR, 9, 1)

  do order = 2, 10    
    YR(order - 1) = REAL(order / 3.e0)
  end do
  call outLines(XR, YR, 9, 2)

  return
end subroutine displayScalingExponent


! ----------------------------
!
! Displays the evolution of all shells in complex space.
!
! ----------------------------
subroutine displayShellValues(VELOCITY, numberOutputs, m)
  implicit none
  
  integer*8 m
  integer*8 n
  integer*8 numberOutputs
  real*4 X(numberOutputs) 
  real*4 Y(numberOutputs)
  complex*16 SHELL(numberOutputs)
  complex*16 VELOCITY(m, numberOutputs)
  
  call PGSUBP(6, 4)
    
  do n = 1, m
    call nShellValues(SHELL, VELOCITY, n, numberOutputs, m)
    call getValuesFromComplex(SHELL, numberOutputs, X, Y)
    call setGraphViewAdjusted(X, Y, numberOutputs, '', '', '', 0)
    call outPoints(X, Y, numberOutputs, 1)
  end do
  
  return
end subroutine displayShellValues


! ----------------------------
!
! Displays the velocity map of the last run.
!
! ----------------------------
subroutine displayVelocityMap(VELOCITY, numberOutputs, m)
  implicit none

  integer*8 m
  integer*8 numberOutputs
  real*4 MODULUS(m, numberOutputs)
  complex*16 VELOCITY(m, numberOutputs)
  
  call shellModulus(MODULUS, VELOCITY, numberOutputs, m)
  call PGENV(0., REAL(m) + 1., 0., REAL(numberOutputs) + 1., 0, 2)
  call outGrayMatrix(MODULUS, m, numberOutputs)
  
  return
end subroutine displayVelocityMap


! ----------------------------
!
! Displays energy through time of the last run.
!
! ----------------------------
subroutine displayEnergy(VELOCITY, numberOutputs, m)
  implicit none
  
  real*8 E
  integer*8 m
  integer*8 n
  integer*8 k
  integer*8 numberOutputs
  real*4 X(numberOutputs)
  real*4 Y(numberOutputs)
  real*8 ENERGY(numberOutputs)
  complex*16 U(-1:m+2)
  complex*16 VELOCITY(m, numberOutputs)
  
  U = (0.d0, 0.d0)
  
  do k = 1, numberOutputs
    do n = 1, m
      U(n) = VELOCITY(n, k)
    end do
    
    ENERGY(k) = E(U, m)
  end do
  
  call getValues(ENERGY, numberOutputs, X, Y)
  call setGraphViewAdjusted(X, Y, numberOutputs, '', '', '', 0)
  call outLines(X, Y, numberOutputs, 1)
  
  return
end subroutine displayEnergy


! ----------------------------
!
! Displays shells modulus through time of the last run.
!
! ----------------------------
subroutine displayShellsModulusThroughTime(VELOCITY, numberOutputs, m)
  implicit none
  
  integer*8 m
  integer*8 n
  integer*8 k
  integer*8 option
  integer*8 numberOutputs
  real*4 X(m)
  real*4 Y(m)
  real*8 MODULUS(m)
  complex*16 VELOCITY(m, numberOutputs)
  character*(10) string
 
  string = ""
     
  print *, 'Type of values:'
  print *, '1  Normal'
  print *, '2  Logarithmic'
  read *, option
  
  k = 1
  do while (k <= numberOutputs .and. string == "")
    call initPGPlot()
    do n = 1, m
      MODULUS(n) = REAL(CDABS(VELOCITY(n, k)))
    end do
    
    if (option == 1) then
      call getValues(MODULUS, m, X, Y)
      call PGENV(0., REAL(m) + 1., 0., 1., 0, 0)
    else if (option == 2) then
      call getLogValues(MODULUS, m, X, Y)
      call PGENV(0., REAL(m) + 1., -30., 1., 0, 0)
    end if
      
    call outLines(X, Y, m, 1)
    read "(a)", string
    call PGEND
    k = k + 1
  end do
  
  return
end subroutine displayShellsModulusThroughTime


! ----------------------------
!
! Displays each shell modulus through time of the last run.
!
! ----------------------------
subroutine displayEachShellModulus(VELOCITY, numberOutputs, m)
  implicit none

  integer*8 k
  integer*8 n
  integer*8 m
  integer*8 starts
  integer*8 length
  integer*8 option
  integer*8 numberOutputs
  real*4 X(numberOutputs)
  real*4 Y(numberOutputs)
  real*8 MODULUS(numberOutputs)
  complex*16 VELOCITY(m, numberOutputs)
  character*(10) string
  
  string = ""
   
  print *, 'Starts at'
  read *, starts
  print *, 'Sequence length'
  read *, length
  print *, 'Type of values:'
  print *, '1  Normal'
  print *, '2  Logarithmic'
  read *, option
  
  do while (starts < numberOutputs - length .and. string == "")
    starts = starts + length / 10
    call initPGPlot()
    call PGSUBP(2, m / 2)
    
    do n = 1, m
      do k = 1, numberOutputs        
        MODULUS(k) = CDABS(VELOCITY(n, k))
      end do
    
      if (option == 1) then
        call getValues(MODULUS, numberOutputs, X, Y)
        call PGENV(REAL(starts), REAL(starts + length), 0.0, .5, 0, 0)
      else if (option == 2) then
        call getLogValues(MODULUS, numberOutputs, X, Y)
        call PGENV(REAL(starts), REAL(starts + length), -20., 0., 0, 0)
      end if
      
      call outLines(X(starts:starts + length), Y(starts:starts + length), &
                    length, 1)
    end do
    
    read "(a)", string
    call PGEND
  end do
  
  return
end subroutine displayEachShellModulus


! ----------------------------
!
! Computes ands prints the correlation of the modulus between each pair of
! shells.
!
! ----------------------------
subroutine shellsCorrelations(VELOCITY, numberOutputs, m)
  implicit none
  
  real*8 Correlation
  real*8 CORRTABLE(m, m)
  integer*8 k
  integer*8 n
  integer*8 m
  integer*8 numberOutputs
  complex*16 VELOCITY(m, numberOutputs)
  
  do n = 1, m
    do k = 1, m
      CORRTABLE(n, k) = Correlation(VELOCITY, n, k, numberOutputs, m)
    end do
  end do
  
  call printMatrix(CORRTABLE, m, m)
  
  return
end subroutine shellsCorrelations


! ----------------------------
!
! Prints the values of the matrix given
!
! ----------------------------
subroutine printMatrix(MATRIX, sizeX, sizeY)
  implicit none
  
  integer*8 i
  integer*8 j
  integer*8 sizeX
  integer*8 sizeY
  real*8 MATRIX(sizeX, sizeY)
  
  do i = 1, sizeY
    write (*, '(a,$)') '['
    do j = 1, sizeX
      write (*, '(a2,f4.2,$)') ' ', MATRIX(j, i)
    end do
    
    print *, ']'
  end do
  
  return
end subroutine printMatrix


! ----------------------------
!
! Returns the correlation of the shell modulud between the shells i and j
!
! ----------------------------
real*8 function Correlation(VELOCITY, i, j, numberOutputs, m)
  implicit none
  
  real*8 si
  real*8 sj
  real*8 covij
  real*8 averagei
  real*8 averagej
  integer*8 k
  integer*8 i
  integer*8 j
  integer*8 m
  integer*8 numberOutputs
  complex*16 VELOCITY(m, numberOutputs)
  
  si = 0.d0
  sj = 0.d0
  covij = 0.d0
  averagei = 0.d0
  averagej = 0.d0
  
  do k = 1, numberOutputs        
    averagei = averagei + CDABS(VELOCITY(i, k))
    averagej = averagej + CDABS(VELOCITY(j, k))
  end do
  
  averagei = averagei / numberOutputs
  averagej = averagej / numberOutputs
  
  do k = 1, numberOutputs        
    si = si + (CDABS(VELOCITY(i, k)) - averagei)**2
    sj = sj + (CDABS(VELOCITY(j, k)) - averagej)**2
  end do
  
  si = SQRT(si / numberOutputs)
  sj = SQRT(sj / numberOutputs)
  
  do k = 1, numberOutputs        
    covij = covij + (CDABS(VELOCITY(i, k)) - averagei) * &
                    (CDABS(VELOCITY(j, k)) - averagej)
  end do
  
  covij = covij / numberOutputs
  Correlation = covij / (si * sj)
  
  return
end function Correlation
